library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(sdmTMB)
library(tidyr)
Mesh
# 300 knots
plot(mesh)
Model specifications
m1 <- sdmTMB(
bio.adult ~ 0 + YearF + GearCat + s(Month, bs = "cc", k = 10),
data = dat,
mesh = mesh,
offset = log(dat$SweptArea),
family = tweedie(link = "log"),
spatial = "on",
time = "Year",
spatiotemporal = "IID"
)
## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.8.1
## Current TMB version is 1.9.0
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
## Warning: The model may not have converged: non-positive-definite Hessian matrix.
Some diagnostics
dat$resids <- residuals(m1) # randomized quantile residuals
hist(dat$resids)
qqnorm(dat$resids[is.finite(dat$resids)])
qqline(dat$resids[is.finite(dat$resids)])
# check spatial autocorrelation
ggplot(filter(dat,Year%%5==1), aes(lon, lat, col = resids)) +
scale_colour_gradient2() +
geom_point() + geom_sf(data=map,inherit.aes = FALSE)+xlim(-12,30)+ylim(47,62)+
facet_wrap(~Year,ncol=2)
## Warning: Removed 758 rows containing missing values (geom_point).
Terms
m1
## Spatiotemporal model fit by ML ['sdmTMB']
## Formula: bio.adult ~ 0 + YearF + GearCat + s(Month, bs = "cc", k = 10)
## Mesh: mesh
## Time column: Year
## Data: dat
## Family: tweedie(link = 'log')
##
## coef.est coef.se
## YearF1967 0.00 NaN
## YearF1968 0.00 NaN
## YearF1970 11.34 NaN
## YearF1971 10.15 NaN
## YearF1972 10.10 NaN
## YearF1973 10.01 NaN
## YearF1974 9.99 NaN
## YearF1975 9.08 NaN
## YearF1976 10.04 NaN
## YearF1977 9.53 NaN
## YearF1978 11.58 NaN
## YearF1979 11.61 NaN
## YearF1980 12.24 NaN
## YearF1981 13.35 NaN
## YearF1982 11.87 NaN
## YearF1983 12.14 NaN
## YearF1984 12.45 NaN
## YearF1985 11.03 NaN
## YearF1986 11.23 NaN
## YearF1987 10.88 NaN
## YearF1988 11.42 NaN
## YearF1989 12.21 NaN
## YearF1990 10.85 NaN
## YearF1991 10.92 NaN
## YearF1992 10.67 NaN
## YearF1993 11.15 NaN
## YearF1994 10.88 NaN
## YearF1995 10.50 NaN
## YearF1996 10.38 NaN
## YearF1997 10.63 NaN
## YearF1998 10.63 NaN
## YearF1999 10.60 NaN
## YearF2000 10.32 NaN
## YearF2001 10.08 NaN
## YearF2002 9.61 NaN
## YearF2003 9.65 NaN
## YearF2004 9.71 NaN
## YearF2005 9.46 NaN
## YearF2006 9.57 NaN
## YearF2007 9.93 NaN
## YearF2008 10.16 NaN
## YearF2009 9.47 NaN
## YearF2010 9.92 NaN
## YearF2011 10.17 NaN
## YearF2012 9.84 NaN
## YearF2013 9.72 NaN
## YearF2014 9.80 NaN
## YearF2015 9.81 NaN
## YearF2016 9.93 NaN
## YearF2017 9.59 NaN
## YearF2018 9.16 NaN
## YearF2019 8.59 NaN
## YearF2020 8.47 NaN
## YearF2021 8.32 NaN
## GearCatBT 0.25 NaN
## GearCatGOV_CL 1.54 NaN
## GearCatGOV_GG 0.97 NaN
## GearCatOTT -3.26 NaN
## GearCatPEL -2.05 NaN
## GearCatTV 0.89 NaN
##
## Smooth terms:
## Std. Dev.
## sds(Month) 0.43
##
## Dispersion parameter: 665.08
## Tweedie p: 1.61
## Matern range: 1.35
## Spatial SD: 0.37
## Spatiotemporal SD: 0.26
## ML criterion at convergence: 248325.051
##
## See ?tidy.sdmTMB to extract these values as a data frame.
#g <- visreg::visreg(fit, xvar = "GearCat", gg = TRUE)
#plot(g)
p1=visreg::visreg(m1, xvar = "GearCat", gg = TRUE)
p2=visreg::visreg(m1, xvar = "YearF",gg = TRUE)
p3=visreg::visreg(m1, xvar = "Month",gg = TRUE)
p1
p2
p3
Predict onto the grid
grid=grid %>% filter(!Year %in% c(1967,1968))
grid$Year=as.integer(grid$Year)
grid$YearF=as.factor(grid$Year)
dat %>% group_by(GearCat) %>% summarise(n=n())
## # A tibble: 7 × 2
## GearCat n
## <fct> <int>
## 1 BAK_GG 520
## 2 BT 33421
## 3 GOV_CL 37009
## 4 GOV_GG 6353
## 5 OTT 87
## 6 PEL 88
## 7 TV 10896
grid$GearCat = 'GOV_CL'
# calculate the mode of month
# getmode=function(x){uniq=unique(x)
# uniq[which.max(tabulate(match(x,uniq)))]}
# getmode(dat$Month) Month set to 9
grid$Month = 9
p <- predict(m1, newdata = grid)
Spatial random fields
plot_map(p, "omega_s") +
scale_fill_gradient2() +
ggtitle("Spatial random effects only")
## Warning: Removed 4817 rows containing missing values (geom_raster).
Spatio-temporal random fields
plot_map(filter(p,Year%%5==1), "epsilon_st") +
scale_fill_gradient2() +
facet_wrap(~Year,ncol = 2) +
ggtitle("Spatiotemporal random effects only")
## Warning: Removed 1169 rows containing missing values (geom_raster).
plot_map(filter(p,Year%%5==1), "exp(est)/1000") +
scale_fill_viridis_c(
trans = "pseudo_log",
# trim extreme high values to make spatial variation more visible
na.value = "yellow", limits = c(0, quantile(exp(p$est)/1000, 0.995)),'cod adult biomass \n (kg km^-2)'
) +
ggtitle("Prediction (fixed effects + all random effects)",
subtitle = paste("maximum estimated biomass density (whole time series) =", round(max(exp(p$est))/1000))
)+facet_wrap(~Year,ncol=2)
## Warning: Removed 1169 rows containing missing values (geom_raster).
